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1. Introduction 

Since its inception in 1992 by White|jl|, the Density Matrix Renormalization Group 
(DMRG) has emerged as one of the most powerful numerical methods in the study of 
low dimensional strongly correlated fermionic, bosonic or quantum magnetic systems 
The universality of its core idea of deducing a decimation prescription for state spaces 
by considering the renormalization flow of suitable reduced density matrices has led 
many researchers to extend the method to other fields of research in systems with many 
correlated degrees of freedom. 

Originally, the method was applied to the renormalization of effective low-energy 
Hamiltonians to study static and dynamic T = properties. Major progress occurred 
with Nishino's realisation^ that the DMRG can be used to renormalize the transfer 
matrix of semi-infinite two-dimensional strips, a method which we will refer to as TMRG 
(transfer matrix renormalization group) in the following. The use of a transfer matrix 
implies that the system is truly infinite in one spatial direction, whereas the power of the 
DMRG to treat large one- dimensional systems is reflected in the fact that the strip width 
can be chosen so large as to allow reliable finite-size extrapolations. In analogy to the 
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approach used in Quantum Monte Carlo, one- dimensional quantum systems were then 
mapped by a Trotter-Suzuki checkerboard decomposition to a two-dimensional classical 
systemQ. The renormalization of the resulting quantum transfer matrix allows the 
study of the thermodynamics of quantum chains for infinite system sizes and very fine 
Trotter-Suzuki decompositions, i.e. down to very low temperatures. 

As the Trotter-Suzuki decomposition generates a transfer matrix which is not 
invariant under a spatial reflection of the system, the quantum transfer matrix is 
asymmetrical and there are left and right eigenstates {ipil and associated with 
the thermodynamically relevant largest eigenvalue of the transfer matrix. This raises 
the question which is the correct way to construct the optimal density matrix for the 
DMRG. In the usual hermitean DMRG, the density matrix is typically obtained by 
tracing out a part of the universe ("environment") in the ground state: 

Psys = Trenv (l) 

In the nonhermitean case, Bursill et al 0] made the symmetrical choice 

Psys = Trcnv (2) 

which can be easily extended to 

P:^r = ^^en.{\^R){U + \^L){^L\) (3) 

which treats the left and right eigenvectors on an equal footing (all composite density 
matrices have to be appropriately normalized, depending on the normalization chosen 
for the eigenvectors; we do not show these factors explicitly in the paper). This approach 
also has the property to minimize the sum of the two squared distances between 
the original eigenstates and their projections onto the reduced state spaces after the 
renormalization step, which is the prescription from which the DMRG procedure for 
symmetrical matrices can be derived. 

However, it was pointed out by Wang and Xiang||^ that the choice of an 
asymmetrical density matrix, 

p:;r = Trenv|V^i?)(V'L|, (4) 

is more physically adequate and yields numerically more satisfying results. This is now 
the accepted choice for the application of the TMRG to quantum problems[^ ^, ||. 

A further field of applications was opened up by the application of the DMRG 
idea to the renormalization of genuinely asymmetrical problems & |10|, O] as they occur 



in nonequilibrium statistical mechanics. Carlon, Henkel and Schollw6ck|ll] exploited 
the fact that the time evolution of react ion- diffusion systems can be mapped to a 
Schrodinger-like equation. In this case, one can study the long-time behaviour of finite 
size systems quite precisely 0, O, IT| . 



Here, the question of the correct choice of the density matrix arises as the transition 
matrix is genuinely asymmetrical itself, since there is no detailed balance. Carlon et al 
found the choice of the symmetrical density matrix to be most suitable, but essentially 
due to reasons of numerical stability: in this class of problems, very high numerical 
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precision for the chosen basis states in the reduced basis was found to be important, 
which is less easily obtained for asymmetrical matrices. An attempt to include projectors 
like \iPr){iPl\ while building a symmetrical density matrix by studying 

p-ir = TrUM + \i^R))im + {^rD, (5) 

yielded inferior results. 

Recently, the TMRG has been used by Kemper et al to successfully renormalize 
stochastic transfer matrices [|15|, an approach which henceforth we will refer to as 
"stochastic TMRG". In their study of the Domany-Kinzel cellular automaton, they 
found that the use of a symmetrical density matrix yielded satisfactory results, while the 
asymmetrical choice led to incorrect results, which is at variance with the conventional 
(or "quantum") TMRG. These findings were obtained on an entirely empirical basis. 

In this paper, we want to investigate in depth the question of the correct choice 
of the density matrix, which is at the core of all asymmetrical DMRG applications, for 
the new stochastic TMRG. We will give explicit constructions for the eigenvalues and 
selected eigenvectors of the stochastic transfer matrices, highlighting the central role of 
the boundary conditions in time direction. This will allow us to demonstrate that the 
content of the asymmetrical density matrix is essentially trivial, while the symmetrical 
density matrix is the physically adequate choice. 



2. Explicit construction of the unrenormalized TMRG transfer matrix 

2.1. The TMRG transfer matrix 

In the stochastic TMRG (for a description of the very similar TMRG applied to quantum 
systems, see 0), one considers a stochastically evolving system extended infinitely in 
one spatial dimension with (for simplicity) local update rules involving neighbouring 
sites only that allow a finite number of states (n). The local interaction between two 
(or similarly, a few) lattice sites is given by the local transfer matrix 

Ut ^r. 



WK = [exp(At.i/i_i)]K?- ' 



fAt 







(6) 



where the stochastic "Hamiltonian" -ffiocai gives the transition rates between states 
and the arrow indicates the direction of the time step At. We assume spatial parity 
invariance 

ir)tl = {r):il . (7) 

In complete analogy to Quantum Monte Carlo and the conventional TMRG, the real 
time evolution operator for the full lattice is mapped approximately by a Trotter- Suzuki 
decomposition to a checkerboard structure of local interactions (figure]^). 

In the time direction the evolution of the system is simulated during a finite 
interval. As illustrated in figure |I|, the full evolution operator is then written as 
exp(2MAt ■ H) = lim^r^oo^^; where M is the Trotter number, N the (diverging) 
system size and H = limN^oo J2iL-N Hiocaiihi + !)■ T is called the transfer matrix 
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Figure 1. Trotter-Suzuki decomposition of the full time evolution into local 
interactions r (2M = 6 time steps shown). The lattice is spatially infinitely extended, 
while in time direction the system is finite. Two adjacent columns (boxed) constitute 
the basic building block: the transfer matrix (T)ri.'.'.r2"M~-i ■ 



which consists of two adjacent columns of r's, where we take as upper (lower) index 
the "history" (time evolution) of the lattice site at the left (right) side of the zig-zagged 
column: 

M-l 

\hk+2^2k+2 



mora 



mi...m2M-i l2M'm.2M k=0 



+ini2fe+i 



(8) 



open initial b.c. open final b.c. 

where /, m, and r denote the left, middle and right lattice sites of T, respectively. One 
uses stochastic initial conditions: the weight of each possible initial configuration is 
given by the product of local weights w{si) of the states Sj at each lattice site i, 



W[. . . ,Si,S2,S3, . . .J 



J|w(si), J2 '^(^i) = 1 (norm.' 



(9) 



Si = l 



such that if each state is weighted equally, there is no bias in the initial configuration. 
At the end, one uses open final boundary conditions, i.e. all final states are allowed 
without bias. In the diagrams, we thus implicitly trace over all internal indices, sum 
over all final indices, and perform a weighted sum over all initial indices. 



2.2. Spectrum and eigenvectors of the transfer matrix 

As the local stochastic transfer matrix r must obey probability conservation, the 
probability of ending up in any state is unity: 



hr2 
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(10) 
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This implies 



TrT 



(11) 



where tracing contracts the upper index of T with the lower one; thus the probability 
that the lattice sites at the left and right boundary of T have the same time evolution 
is unity. In pictorial language, the trace rolls up the transfer matrix into a cylinder, 
with its left and right boundary identified and the middle site at the opposite side of 
the cylinder (cf. figure H). 
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Figure 2. Tr{T)r\::.?2M^i is equivalent to T being rolled up on a cylinder with a solely 
local interaction. 



But this is equivalent to having only the local interaction between the left/right 
lattice site on one side and the middle one on the other, for 2M time steps At, which is 
known exactly (e.g. 2M = 6): 
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(12) 



moro leme h...l5mi...ni5 

r^^'^(At) (by parity invar. (7)) 

J r'omo 

= Y: ^(mo) w{ro) E [exp(2MAt ■ i/iocai)]':ro 

moro iQmQ 
moro iQiriQ 



Y w{mo) w{ro) -1 = 1 (by (|T0|) and 



moro 



Next we consider the spectrum and eigenvector decomposition of T. As all entries 
of T are non-negative, the Frobenius theorem implies that the largest eigenvalue Aq 



Density matrix in the stochastic TMRG 



6 



is non-degenerate and real. The associated left and right eigenvectors {i/jlI and lipn) 
(where \iPr)^ 7^ {ipil in general, since T is not symmetrical) have nonnegative entries. 
Now the decisive difference between conventional and stochastic TMRG is that in the 
conventional TMRG one has periodic boundary conditions in (imaginary) time, but 
open final boundary conditions in the stochastic TMRG. This latter property implies 

T^*^-! = \iIjr){iPl\ = T^ VA; > 2M - 1 , (13) 

i.e. sufficiently large powers of T are identical and decompose into an outer product of 
the two eigenvectors of Aq = 1, normalized as {ipLl'ipR) = 1, while all other eigenvalues 
of T vanish. 




i 4 



Figure 3. Open b.c. lead to trivial t's that cause also the final indices of the previous 
t's to be summed over, thus "propagating" the open b.c. backwards in time. 



To see this, consider figure ^ Because r is stochastic, it becomes trivial by summing 
over its final indices (cf. (p!OD; its value does not depend on its initial indices, such that 
the contraction with the final index of a previous t yields an unweighted sum over the 
previous r's final index. Therefore in figure ^, the trivial r's "propagate" the open final 
b.c. backwards in time, resulting in a cone of trivial r's that have no influence on the 
value of the transfer matrix T''. 
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Figure 4. decompose into a left and right part of ever smaller rank until 

rp2M-i decomposes into two lightcones which cannot interact during 2M time steps 
(here 2M = 4): (a) T = • (V'l)"'^™! ^i^^ 2 connecting indices (mi, ma), 

(b) T2 = {iPr)I^Y'' ■ {■>pL)T,r2r3 '^1*^ ^nc Connecting index (mi), (c) T^^'^-i = 
(^^j^yihh . (V'i)rj^2r3 Hr){iPl\ with no connecting index. 



The Trotter-Suzuki decomposition introduces a finite speed of action ("speed of 
light") into the model: during each time step, only nearest neighbours that are connected 
by an edge of a [nontrivial) t can interact, so interactions can only propagate at a speed 
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of one lattice site per time step At. Hence we speak of the "light cone" of a certain site 
as those sites in time and space that influence the probabilities of different occupations 
for this site. 

In T (cf. figure ^a), the last r is trivial, hence it has the value 1 indepedent of I3 and 
777,3. Thus we effectively sum over the internal index 777.3, and ^3 is a "dummy index": 
the value of iT)l}^^^l^ does not depend on it. This leaves only two indices (7771,7772) 
connecting the left and right part of T, so the transfer matrix can be decomposed into 
a product of two matrices ipR and ipL of rank 77^ (77 is the number of possible states per 
site) . 

When we multiply this with another T (cf. figure ^o), we observe that a third r 
becomes trivial as the final r's propagate the final open b.c. backwards. This leads to 
a sum also over 7772, leaving 7771 as the only connecting index, so decomposes into a 
product of two matrices of rank 77. 

Finally in (cf. figure ^), the trivial r's propagate all the way to the initial r's, 

effectively separating the left and right boundary: now there are no connecting indices 
left between ipn and ipi, turning them into column and row vectors lippi) and {iPlI- From 
the viewpoint of causality, the left and right boundary are separated by 4M — 2 lattice 
sites, so no initial lattice site can influence both the left and the right boundary, i.e. the 
light cone for the lattice sites on the left boundary (/i, . . . , /2M-1) does not intersect the 
light cone for the lattice sites on the right boundary (ri, . . . ,r2A/-i)- This means that 
the probability of a particular history on the left boundary is completely independent 
of the probability of a particular history on the right. 

We can now explicitly compute these column and row vectors l'?/;^)'^ - '^*^-! and 
(V'i|ri...r2ii/_i- we take the left hght cone (of triangular shape, delimited by the dashed 
line) and contract the included nontrivial local transfer matrices r in all internal, initial, 
and final indices, keeping the free indices li . . . I2M-1, and denote this object as {right 
eigenvector, see below) lipnY^-'-^'^'''^ . Analogously, we denote the same construction on 
the right with free indices ri . . .r2M-i as {4'L\ri...r2M-i- Note that the last index on 
the left, I2M-I1 is by construction a dummy index which does not change the value of 

If we consider even higher powers of T than 2M — 1, we only add more trivial r's 
in the middle that lie outside either light cone and therefore do not change the value of 
the transfer matrix. Thus we have proven T"^^^'^ = \/k > 2M — 1, which implies for 
the eigenvalues Aj of T that A,^^'^^^ = A^*^ Vi, such that Aj can only assume the values 



or 1. From (y) follows Aq = 1, Aj = Vz > 0. 

The above decomposition of into an outer product of two vectors \iPr){4'l\ 

is in fact the explicit eigenvector decomposition corresponding to the eigenvalue Aq = 1: 
multiplying the right (left) eigenvector onto from the right (left) is, in pictorial 

language, attaching the respective light cone to the transfer matrix with a trace over 
all overlapping indices. But now the final open b.c. can propagate through all r's on 
that side, turning them trivial (which automatically gives the correct normalization 
(V'lI'^j?) = 1) and leaving just the light cone on the other side, which reproduces 
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exactly the previously attached eigenvector. The eigenvalue 1 is determined by this 
normalization condition: 

T2^-VH) = |^ii)(^L|^ii) = |^ii)-l. (14) 

This concludes the proof of ([T3|). Similarly, we check that lipR) and are also the 
eigenvectors of T corresponding to the same eigenvalue (cf. figure 
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Figure 5. \iPr) is indeed the right eigenvector of T corresponding to the eigenvalue 

Ao - 1: for 2M - 4, Tll^,t\^,\^ Rr'""'""' = \^rY''''''- 



Let us finally emphasize that for small powers of T, with k < 2M — 1, this 
simple decomposition does not hold. But this is not relevant, as we take k to infinity in 
the TMRG. 

To arrive at a prescription for the stochastic TMRG density matrix, we must 
now consider how actual physical information is extracted from the stochastic transfer 
matrix. This is done by computing expectation values of local operators, averaged over 
the whole lattice at final time 2M ■ At by multiplying r in the last time step with the 
local operator, e.g. for the local particle number operator n{i) one has 

irn)tl := ir)\i:i . nik) (15) 

where we denote the transfer matrix with its last r replaced by r„ as T„ (coloured black 
in the diagram). Then the expectation value of the particle number operator can be 
calculated as 



Tr(r^T„) 
{nt=2M-At) = lim ^ /^Af ,?x = lim 

N^oo Tr(T^^+^) N-*oo 





h d 




Tr(|^«)(^L|) 




(16) 



(17) 



where we have omitted local transfer matrices that become trivial by summing and used 
the pictorial construction of the eigenvectors. The computation of expectation values 
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Figure 6. The light cones that are needed to compute the expectation value of the 
particle number operator for the first few time steps: (a) 2M = 2, (b) 2M — 3, (c) 
2M = 4. 



is thus reduced to the computation of the hght cone of lattice sites that can possibly 
influence the particular site at which the operator is measured. 

Figure |^ shows the light cones for the first few time steps. In this manner one can 
compute the exact evolution of {n{t)), the only inaccuracy stemming from the Trotter- 
Suzuki decomposition. It is interesting to observe that for determining both eigenvectors 
and expectation values, a matrix diagonalization can be completely avoided. 

For the choice of the density matrix the important observation in figure ^ is now 
that by parity invariance (|^ and adding one dummy index {i2M-i), 

l^i?)2^""-^ = ((^l|2M-m,....,,,,„,,.J^ =: \i^L)^M-r' (18) 

i.e. the left eigenvector for a given time step is identical to the right eigenvector in the 
next (half) time step. 

3. Density matrix renormalization of the stochastic transfer matrix 

We have seen how, in principle, physical information can be extracted from the stochastic 
transfer matrix, but both memory usage and computation time grow exponentially 
with the number of time steps. The transfer matrix operates on the space of possible 
time evolutions of particular sites, which has to be decimated to remain numerically 
manageable. This is the idea of the TMRG. As the probability of a partial time 
evolution is determined not only by its past but also by its future continuation, and 
the infinite future is not yet known, one contends oneself with a "lookahead" of several 
time steps which is called the "environment" , while the past time evolution — that is 
being projected onto a smaller Hilbert space — is called the "system" . 

Traditionally, one obtains this basis by computing the reduced density matrix for 
the system, psys, by taking a partial trace over the environment of the outer product of 
right and left T eigenvectors; the new Hilbert space basis consists of Psys's eigenvectors 
associated to the largest eigenvalues (highest probabilities). Inspired from traditional 
TMRG (with periodic b.c. in time direction), one could try Psys™ = Trgnv |V'-r)(^l|- But 
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in the case of open b.c, this is just a partial trace of the transfer matrix, 



asym 
r sys 



Trenv \i^R){i^L 



Tr T 



2M-1 



|V^i?)sys(^ 



L|sys 



(19) 



having a spectrum of 



.3,, , {1,0,..., 0}. (20) 

We demonstrate this for 2M = 4, taking e.g. the first index as the system index and 
the following two indices as environment indices. Then the trace over the environment 
Trenv = J2i3r3i2r2 ^h^h^ iUustratcd in figure ^ makes all r's in the environment trivial: 
when we contract the dummy index Z3 with r^, the final r of {ipil becomes trivial, so r2 
now becomes a dummy index. This in turn gives a sum over I2, making also the final 
r of lipn) trivial, etc. In the end we are left with {iPr) and {iJjlI of an earlier time step, 
namely the one with only the system indices /i, ri. 
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Figure 7. Computing tlie asymmetrical density matrix pfy^™: the trace makes the 
environment trivial, Trenv |V'-r)(V'l| = |''/'-R)sys(V'L|sys. 



Thus this asymmetrical density matrix provides only a single basis vector for the 
Hilbert space projection, which is obviously not enough. The reason is that the open b.c. 
"average away" the physical interaction in the environment by making it irrelevant how 
the system will react in the future. In these boundary conditions resides the essential 
difference between quantum and stochastic TMRG. 

Each renormalization step using Pgyg™ projects therefore onto a one- dimensional 
Hilbert space. No longer-ranged correlation implicit in differently weighted basis vectors 
is being built up. We therefore expect to observe only the local evolution of the model, 
which is exactly what we got from numerical simulation using the asymmetrical density 
matrix. 

As an example, we use an exactly solvable reaction-diffusion model: AA 00 
with rate 2a (reaction) and AQ ^ OA with rate D (diffusion). For 2a = D, the exact 
solution is given by Spouge|T6|: for initial partical number n{t 

w{0) 



0) 



(unbiased 



w{A) 



2' 



i.e. all lattice sites have an independent probability ^ of occupancy). 



(n(t)) = iexp( 



-2Dt) [Io{2Dt) + Ii{2Dt)] ~ ^ (asympt.) 



(21) 



This agrees well only with the simulation using the symmetrized density matrix p^ymm 



sys 



(cf. figure ||). But when we instead use the asymmetrical density matrix p^ys™; 
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observe a 1/t decay. This resembles very closely what we expect from a local interaction: 
the exact 4-site interaction (i.e. renormalizing as soon as the system has n states) is 

{n{2M + 1)) = {n{2M)) - [1 - exp(-DAt)] ■ {n{2M)f (22) 

(n(0) = |) which goes asymptotically as 1/t. 




time t=2MAt 

Figure 8. Comparing the simulation of a reaction-diffusion model for symme- 
trized vs. asymmetrical density matrices. ( ) exact solution, ( ) exact local 

(4-site) interaction, (O ) simulation, symmetrized density matrix, (<)) simulation, 
asymmetrical density matrix (D = 0.1 in inverse time units, and At =1). 

The origin of this shortcoming becomes clear in the following consideration: figure ^ 
shows the left eigenvectors during consecutive time steps. Suppose we want to find a 
good basis for renormalizing {iPl\2M=5, taking its first two indices (Zi, I2) as system and 
the latter two (Z3, U) as environment. After somehow tracing over the environment, the 
new basis vectors will have indices h, I2, and they will roughly resemble {ipiU, having 
a strong correlation between li and I2 because they are connected by an edge of r. 

Now in order to advance one time step from (V'lIs to O'lIg, we see from figure ^ that 
we need to multiply (V'lIs by the matrix represented by the column Cq. This reveals 
the problem: the first two indices of Cq are not connected by an edge of r and therefore 
completely independent. Thus the renormalization basis obtained from {iPlU with a 
strong correlation between the first pair of indices fits very poorly for renormalizing 
Cq. The states in the basis are essentially orthogonal to those that would describe a 
renormalized Cq well. Looking ahead, C7 which would bring us to {ipili has again a 
structure similar to (V'lIs and could be renormalized well. 

This problem would be remedied by adding to the renormalization basis also basis 
vectors of a form similar to {iPlIa, (tracing over the last index as environment) which 
shows little correlation between the first two indices. Noting that {iPlU = (|V'-R.)5)"^ 
(p!8|) we have to mix left and right eigenvectors, {ipiU and {ipnl^, in order to obtain 
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Figure 9. The evolution of the largest left T eigenvector {i}}l\2M for the first few 
Trotter steps. 



a good basis. Based on the original argument by White about the optimal density 
matrix in DMRG, Carlon et al [^] have shown that in order to express two different 
vectors IV'r) most accurately (in the same basis), one has to use an eigenbasis of 

Psys = Trenv(|'?/'R)('0i?l + I'^/'i) (V'i I ) • The Symmetrical form of this density matrix arises 
because they used as a measure of distance between exact and approximated vectors 
the symmetrical norm IIIV'r) — l"^/?)!!^ + III'^l) " With this same argument, we 

now conclude that this symmetrized density matrix is also optimal in the case of explicit 
eigenvectors. 

This explains why a symmetrized density matrix, e.g. Pgy™™ = Trenv (| (^_r| + 
as used by Kemper et al |]15[, works so well. Indeed, within the precision of 
the method, the exact result for the example considered is reproduced when we use this 
symmetrized reduced density matrix Pgy™™. 

If we consider for our example the eigenvalue spectrum of the symmetrized density 
matrix, initially it falls off approximately exponentially, Aj ~ a~* for some constant a (cf. 
figure p!0| ). After many renormalization steps, the spectrum falls off rather according 
to a power law with large negative exponent, Aj ~ i~^ . The symmetrized density 
matrix therefore indeed fulfills two basic requirements: It should have as many nonzero 
eigenvalues as possible in order to keep rich nontrivial information about the future 
evolution in the environment (contrary to (pp])), but at the same time they should drop 
off as quickly as possible so that when we truncate the spectrum at m states that are 
retained, we lose as little information on the system as possible. 

Now, in practical applications of the quantum TMRG, system and environment are 
treated on an equal basis, renormalizing the system with a density matrix traced over 
the environment and vice versa. In the stochastic TMRG, the explicit construction of T 
and the eigenvectors opens the possibility to renormalize only the system and keep the 
environment as a "lookahead" of constant finite length. This simplifies the algorithm. 
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Figure 10. Spectrum of the reduced system density matrix Pg^™™ at specific steps 2M 
(to = 64). (O ) 2M = 18 (trunc. 10-^^°), {<}) 2M = 20 (trunc. lO^^S), (y) 2M = 60 
(trunc. f0^2i)^ (□) 2M = 280 (trunc. lO^^i). 



since the ^/''s can be renormalized explicitly as shown above, and there is no need to 
ever solve the eigenvalue problem of the transfer matrix; on the other hand, due to 
the exponential growth of the matrices with a longer lookahead, this approach is in 
practice limited (see our numerical example below), and one reverts to the conventional 
approach of renormalizing both system and environment of the transfer matrix T and 
then finding its largest eigenvectors anew in each step by using a suitable algorithm for 
the diagonalisation of large nonsymmetrical matrices. 

Our considerations above were all for the case without renormalized environment. 
Renormalizing it, too, leads to an increasingly longer, but less precise lookahead. 
As our considerations hold formally for all lengths of lookaheads, the result on the 
choice of the density matrix should hold also in that case, which is indeed what we 
observe numerically. To check whether our special choice of a symmetrical density 
matrix is indeed optimal among symmetrical density matrices, we have also considered 
numerically several other symmetrical density matrices for comparison, in particular 

Trenv(|V'L)(^fi| + IV^fl)(^L|), (23) 

Tt,U\^l) + \^R)mL\ + {M). (24) 

Trenv(IV'L)- W)((V^L|-(V'i?l), (25) 

but the convergence turned out to be worse than for Psy™"^. AH results shown in the 
following are therefore for the "canonical" choice of the symmetrical density matrix. 

As an application, we tried the branching- fusing process explained in |Tl|: diffusion 
AO ^ OA (rate D), and several reactions AA 00 (rate 2a), AA OA, AO (rate 7), 
OA, AO —>■ 00 (rate 6), and OA, AO — *• AA (rate (3). The relations between the rates are 
fixed a.sD = 2a = 'j = 6 = l— p, P = p. There is a critical point of the directed 
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percolation universality class at Pc ~ 0.8403578. First we check the dependence of our 
TMRG simulation (using the conventional algorithm) for different values of m (number 



of states kept during projection into a smaller Hilbert space) (cf. figure 0) at the critical 
point, the most difficult point for the stochastic TMRG. 
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Figure 11. Branching-fusing model|Q at the critical point p = 0.8403578 for different 
number of states kept (m) resp. truncation in the spectrum of psys- 



We observe that the particle number generally decays exponentially, despite 
criticality, because we are operating at very short times and due to the rather rough 
Trotter decomposition. On the lower branch, for a fixed Trotter time {At = 1), our 
results become more accurate as we increase m, effectively truncating the spectrum of 
Psys at a smaller eigenvalue, until we reach the machine precision (53 mantissa bits). To 
compute more time steps, one option is to use more accurate floating point numbers 
(see below). 

For comparison, choosing a smaller time step At (i.e. a finer Trotter decomposition, 
last two curves in figure |1T]), while maintaining a large Hilbert space (m = 32), we 
get very different results: the particle number decreases much more slowly, which 
shows that our Trotter decomposition is still far from being sufficiently fine (of course, 
extrapolations At ^ are feasible). Unfortunately the number of time steps we 
can simulate decreases, too, which makes such an extrapolation more difficult. This 
can be explained as follows: the shorter the time step At, the smaller the transition 
probabilities to different states, and the closer is r to the identity matrix. But a transfer 
matrix composed of such r's has eigenvectors with exponentially growing Euclidean 



norm {{^r\^r) 



~ n 



M 



and analogously for This is not changed by renormalization 



because an orthogonal transformation conserves the norm, and we only lose those Hilbert 
space directions that are very improbable. As the renormalization keeps the number 
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of components of {iPl\ and \iPr) bounded, the values of the components have to grow 
exponentially to obtain a growing norm. At the same time, (V'l|^_r) = 1, so {ip^l and 
are almost orthogonal. Now the numerical problem is the following: while the 
components grow exponentially, they have to cancel each other very accurately in order 
to yield the correct normalization. This works only until the components are of order 
of magnitude where e is the machine accuracy. Choosing more accurate floating 

point numbers, by taking twice as many mantissa bits we can simulate very roughly twice 
as many time steps. This is illustrated in figure |ll]: if we use the GNU multiprecision 
library (GMP), truncating at 10~^^ instead of 10"^^, we can simulate approximately 1.5 
times as many steps. 

Finally we want to compare the numerical stability of the conventional TMRG 
algorithm, renormalizing both system and environment, with the newly proposed 
method of explicitly constructing the eigenvectors (saving the iterations of the transfer 
matrix diagonalisation which take up more than half of the time in the original 
algorithm). Here we renormalize only the system, keeping a lookahead of finite length as 
environment. As can be seen in figure |12|, the result is much better than the conventional 
algorithm for an intermediate accuracy of m = 11 states kept, while it cannot compete 
with the most accurate (up to machine precision) values for m = 28. So, as already 
mentioned above, it is more proof of a concept and a tool for analytical derivations than 
for practical simulations. Neither more sites lookahead nor larger m improve the result 
significantly. The computation time for the displayed curves was in the range of a few 
seconds (m = 11) up to a few tens of seconds (m = 32) on a PHI (500 MHz), while 
memory usage was at most a few MBytes. 




Figure 12. Branching-fusing model|0 at the critical point p — 0.8403578 for different 

TMRG methods. ( ) conventional method, m = 28; ( ) conventional method, 

m = 11; (O ) explicit eigenvectors, m = 8 (4 sites lookahead); {(}) explicit eigenvectors, 
TO = 11 (6 sites lookahead). 
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4. Conclusions and Outlook 



We have shown a new way of looking at stochastic TMRG with open b.c. by explicitly 
giving the transfer matrix and its left and right eigenvectors corresponding to the 
physically relevant eigenvalue 1. This allowed us to study the physical effect of various 
choices for the density matrix of the TMRG. The conclusion was that, unlike the TMRG 
for quantum transfer matrices, an unsymmetrical density matrix is not just less accurate, 
but fundamentally unsuitable for the stochastic TMRG. The core problem was that one 
has open boundary conditions in the future of the time evolution, averaging away the 
physical interaction, which is not the case in the periodically closed quantum transfer 
matrix DMRG. The symmetrical density matrix defined by demanding an optimal 
representation of both left and right eigenvectors in the truncated basis turned out 
to be the mandated choice. An important observation was that in certain problem 
classes where small transition probabilities lead to local transfer matrices close to unity, 
there is an inherent problem with the finite precision of computer arithmetic if one 
wants to have fine Trotter decompositions and to access comparatively large real times 
simultaneously. 
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